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We derive a quantum theory of evaporative cooling for a degenerate Fermi gas with two constituents 
and show that the optimum cooling trajectory is influenced signiflcantly by the quantum statistics 
of the particles. The cooling efficiency is reduced at low temperatures due to Pauli blocking of 
available final states in each binary collision event. We compare the theoretical optimum trajectory 
with experimental data on cooling a quantum degenerate cloud of potassium-40, and show that 
temperatures as low as 0.3 times the Fermi temperature can now be achieved. 
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The recent demonstrations of Bose- Einstein condensa- 
tion in dilute alkali and hydrogen gases have required the 
ability to reach extremely low temperatures in the micro- 
Kelvin to nano-Kelvin scale. To date, this has been pos- 
sible only by using the experimental technique of forced 
evaporative cooling . Efficient evaporative cooling can 
allow the temperature of a gas to be reduced by orders 
of magnitude without prohibitive loss in the number of 
atoms. It has universal application to cool magnetically 
trapped atomic and molecular vapors and has already 
been applied to produce quantum degenerate clouds of 
rubidium |2|, sodium 1^, lithium Q], potassium and 
hydrogen |6|. 

For a bosonic gas, cooling can be continued to the 
point where no discernible normal component of the gas 
is present, closely approximating a zero temperature sys- 
tem. Demonstrating the ability to reach this regime has 
been a prerequisite to many of the recent experiments 
on collective effects in these systems. Collective phe- 
nomena that have now been observed include linear re- 
sponse , surface modes Q , and topological excita- 
tions such as vortices [|2]jl^ . A current goal is to observe 
the conjugate low-temperature phenomena in a fermionic 
gas when it is cooled well below the onset of quantum de- 
generacy 

In evaporative cooling, a "cut" is made at a prescribed 
energy and all atoms with energies greater than the cut 
are removed from the system. The remaining atoms will 
rethermalize by collisions to form an equilibrium with 
a lower temperature. The crucial parameter that de- 
termines the timescale for cooling is therefore the rate 
of rethermalization F. For a dilute gas at temperatures 
where quantum statistics do not play a role, rethermal- 
ization is determined by the elastic collision rate, given 
by F = nav, where n is the spatially averaged density- 
weighted density, a is the collision cross section, and v is 
the root-mean-square velocity of the colliding species. In 
a harmonic trap, F may increase as the gas is cooled de- 
spite the obvious reduction in average velocity. The rea- 
son for this is simply that as the cloud cools, the atoms 
fall to the bottom of the trap and become more tightly 
confined, increasing the number density n and more than 
compensating for the loss in energy per particle. Achiev- 



ing this regime, known as runaway evaporation, is typi- 
cally an experimental prerequisite for following a cooling 
trajectory that leads to a quantum degenerate gas. 

For bosons, once the temperature is reduced below the 
critical temperature for Bose-Einstein condensation to 
occur, effects due to quantum statistics assist the evap- 
oration. Consider a typical collision event involving two 
atoms from the normal thermal gas that initially have 
approximately the mean energy of the distribution. The 
presence of a Bose-Einstein condensate modifies the scat- 
tering probability into each possible final state and en- 
hances the likelihood of stimulated scattering of one of 
the atoms into the condensate. The other atom then 
obtains the total energy of the initial pair and can be 
removed by the evaporative cut. Clearly this type of col- 
lision leads to very efficient evaporative cooling. 

The opposite situation is true for fermions. As the 
temperature falls below the Fermi temperature, efficient 
collisions turn off due to Pauli blocking since the states of 
lowest energy become occupied with high probability . 
In this paper we study this effect on the achievable op- 
timum evaporation trajectory. Our calculations are mo- 
tivated by the first application of evaporative cooling to 
produce a quantum degenerate Fermi gas Q . This recent 
experiment has opened the door to the study of Fermi 
statistics in an extremely dilute regime — perhaps even- 
tually allowing for the possibility of investigating Cooper 
pairing and the BCS phase transition in these dilute sys- 
tems. 

At typical temperatures of interest, collisions between 
atom pairs are purely s-wave since the characteristic col- 
lision energies are well below the centrifugal barrier as- 
sociated with channels of non-zero orbital angular mo- 
menta p9| ]. Since for fermions the total wavefunction 
must be antisymmetric with respect to exchange of any 
pair of atoms, s-wave collisions are only possible if at 
least two internal atomic hyperfine states are simultane- 
ously present, or alternatively if sympathetic cooling is 
performed with a distinguishable species, such as a dif- 
ferent isotope or a different element j2^j2l[]. Here, we 
consider the first of these possibilities — a two-component 
Fermi gas. 

The Hamiltonian for this system may be separated into 
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two parts, H = Hq + Hi where Hq is the usual single par- 
ticle energy of the system and Hi describes binary colli- 
sions. The Hamiltonian Hq for a two-component mixture 
confined in a three-dimensional harmonic oscillator is 



Hq — ^ En ( 



at 



(1) 



where we have assumed the potential is identical for both 
species. The summation is taken over the three inte- 
ger components of n = {nx,ny,nz). If the harmonic 
potential is isotropic with oscillation frequency to then 
En = huj{nx + Uy + Tiz + 3/2). The annihilation opera- 
tors for the two components, a„ and 6„, obey the usual 
Fermi commutation relations. Binary collisions are de- 
scribed by the interaction Hamiltonian 



n,m,q^p 

where the matrix element is 
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The oscillator eigenf unctions a„ and Pn form a complete 
orthonormal basis that spans the two-component Hilbert 
space. In calculating the matrix element, we have re- 
placed the physical two-particle potential by a contact 
potential. The dimensional prefactor is Uq — ATrh^a/m, 
where m is the atomic mass, and a is the s-wave scatter- 
ing length, which includes contributions from both direct 
and exchange scattering. 

Although, in principle, one could solve the evolution of 
this isolated many-body system, we are primarily inter- 
ested here in a simplified description on a coarse-grained 
timescale. Such a description is given by quantum kinetic 
theory where a set of relevant observables are quantities 
of interest. The underlying theoretical framework derives 
from the property that collisions in the dilute gas are ex- 
tremely well separated in time. This allows the Born and 
Markov approximations to be made and the subsequent 
derivation of a pcrturbative theory in lowest orders of the 
interaction Hamiltonian ||2^]. The relevant observables 
are the populations of the two species (diagonal elements 
of the single particle density matrix): An = (aJiOri) 

Br. = {bib J. 

This approach must be extended in order to treat 
boson-fermion mixtures, or Fermi gases at temperatures 
where Cooper pairing is important. In those situations it 
would be necessary to expand the set of relevant observ- 
ables to consider diagonal and off-diagonal contributions 
to the normal and anomalous densities, as well as the role 
of mean-fields ||2^ . 

Following this procedure, the quantum kinetic equa- 
tions (also referred to as the "Quantum Boltzmann Equa- 
tions" ) for the two-component Fermi gas are given by 



dAn 

dt 
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where the transition rates are found from Fermi's golden 
rule 
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Note the factors such as (1 — An) and (1 — -Brn) give rise 
to the mechanism known as Pauli blocking. The Kro- 
necker delta function, 5, constrains the collision to be 
precisely on the energy shell. This is an unphysical ar- 
tifact of the Markov approximation that is valid for the 
dilute gas only when the elastic collision rate is much less 
than the oscillation frequency in the trap. Otherwise the 
collisional broadening of the levels would be greater than 
their spacing, and the energy basis we use here would 
not be an appropriate choice since off-diagonal elements 
would then be important. We assume ergodicity by as- 
signing equal population to each state in the degenerate 
manifold of states with the same principal quantum num- 
ber n (and therefore the same energy). We denote the 
ergodic populations by A^^ and B^^ (indexed only by 
the discrete values of the energy e„) which are related to 
previously defined populations for an arbitrary quantum 
state An and Bn by 



n 
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with (7e„ denoting the degeneracy of states for the three- 
dimensional harmonic oscillator 
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(n + l)(n + 2). 
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The quantum kinetic equations given in Eq. (^ may be 
simplified by approximating the summations by integrals 
over continuous distributions. This is typically always a 
good approximation for fermions, since for a sufficiently 
large sample, a macroscopic number of states are occu- 
pied even at very low temperature. The quantum kinetic 
equations then describe the rate of transfer of popula- 
tion between continuous distribution functions. We de- 
note these functions (of a continuous energy variable e) 
as A{e) and B{e), which evolve according to 
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deep(e)F(e) = E. 



(10) 



at TT^Tr J 

^A{ep)B{eq)\l — i?(em)][l — A(e„)] ^j^g truncated Fermi-Dirac distribution function is illus- 

, , , , , ^,1 , ^ trated in Fie. nl for the initial conditions of the evapora- 

-A{e^)B{em)[l - i?(e,)][l - ^(e,)]}, (8) ^i^^i^tion^ 



where emin = min{e„, e^, Cp, e,} and A = e„ + e™ - 
Gq — &p- The quantum mechanical cross-section applicable 
here is cr = 47ra^ since the products of a collision event are 
in quantum mechanically distinguishable spin states [ p4| . 
The density of states p{e) for the three-dimensional har- 
monic oscillator can be found from the large n limit of 
Eq. which gives p(e) = ^e^. Although these equa- 
tions are similar in form to the discrete version in Eq. (Q), 
a key simplification has been the replacement of the col- 
lision kernel [defined in Eq. (H)] by p(einin); which is the 
classical limit As was shown in Ref. the con- 

vergence to the classical limit is very rapid as Cmin is 
raised. Significant quantum correction occurs only when 
both of the colliding atoms are in the lowest few (of or- 
der one to five) states of the harmonic oscillator. These 
lowest energy collisions give a microscopic correction to 
the collision rate when the particles have Fermi statistics 
and are typically distributed over a macroscopic number 
of levels of the oscillator. The situation is quite differ- 
ent for the quantum degenerate Bose gas where a careful 
treatment of these low energy collisions is usually crucial 
due to the possibility of condensate mean-fields. 

While the simultaneous equations in Eq. (||) may be 
solved by direct numerical integration, the calculation 
is cumbersome given the multidimensional integrals that 
must be performed at each timestep. We provide a more 
intuitive approach that is motivated by the near equi- 
librium distributions expected when the elastic collision 
rate is sufficiently high. We assume that the form of the 
population distribution functions for both components 
are given by a truncated Fermi-Dirac distribution F{e) 
as defined by 
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e<K 
otherwise 
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A similar method was introduced in Ref. ||2^ to treat 
the evaporative cooling of a classical gas. For simplicity, 
we also take the simplest case of ^(e) = B{e) = F{e) 
at all times, since if the distributions of the two com- 
ponents are initially identical they will remain so due 
to the symmetry of the equations. This means that the 
distribution functions for both species are parameterized 
by the same three variables: (i) an inverse temperature 
/?, (ii) a chemical potential /i, and (iii) a cut energy K. 
Given a value for the cut energy, (3 and p can be solved 
to give simultaneously the correct total number of atoms 
in one component N , and total energy of these atoms E, 
according to 
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FIG. 1. The truncated Fermi-Dirac distribution. The 
graph shows the distribution function, p{e)F{e), as a function 
of e/fcfl where fcs is Boltzmann's constant. The parameters 
used were those for the start of the evaporation simulation 
with iV = 3.2 X 10^ E = SNksT where T = 1.5 ^tK, and 
K — 4E/N. The inset shows a magnified view of the discon- 
tinuity at the cut energy. 

The simulation algorithm is the following: 

1. Starting with given N, E, and K, solve Eq. ([l^) to 
find /? and p. 

2. Consider propagation of the kinetic equations for a 
time step dt and determine the change in number 
dNi and energy dEi due to atoms colliding and 
gaining energy above the cut K , 



dNi = dt 
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dEi = dt 
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3. Simulate background loss (energy-independent re- 
moval of atoms arising from non-ideal vacuum con- 
ditions in experiments) with rate 7 from the trap 



dN2 'jNdt, 
dE2 = jEdt. 



(12) 



4. Lower the cut energy, from to a new value K' 
and find the change in number and energy due to 
trimming the highest energy atoms: 



3 



dN3= [ dep(e)F{e), 

JK' 

dEs = [ deep{e)F{e). 

JK' 



(13) 



5. Update the number and energy 
N ^ N -^dN^, 

(T 

E ^ E -^dE^, 



(14) 



and repeat this sequence starting again from step 1. 

A technical point is that solving Eq. (^^ for f3 and 
/i is a two-dimensional root finding problem that can 
potentially be non-trivial. We use a multidimensional 
Newton-Raphson algorithm that will rapidly produce a 
good estimate of the value of the solutions in a few iter- 
ations. To find a good estimate to start with, we employ 
a simple three-point polynomial extrapolation of solu- 
tions from previous timesteps. In this extrapolation we 
use the cut energy K as the independent variable, rather 
than the step number or time. 

Although this method may be used to calculate the 
evaporation trajectory for an arbitrary time dependence 
of the cut energy, we are most interested here in deter- 
mining the optimum path. During the evaporation simu- 
lation, we follow a trajectory that maximizes the energy 
removed per particle from the system. That is, we choose 
a value for K' in such a way as to numerically maximize 



(15) 



for the subsequent timestep. 

In Fig. 1^, we show the calculated optimum evaporation 
trajectory for the two-component Fermi gas. In order to 
indicate the level of quantum degeneracy, we have nor- 
malized energies and temperatures by dividing them by 
the Fermi energy and Fermi temperature (see caption). 
The optimum cut energy approaches closely the Fermi 
surface towards the end of the simulation. While the 
ideal evaporative trajectory demonstrates the theoretical 
possibility for achieving very low temperatures — with the 
chemical potential tending towards the Fermi energy and 
with a macroscopic population remaining — the efficiency 
of the evaporation trajectory falls dramatically as the 
system becomes degenerate. 
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FIG. 2. Optimum evaporation trajectory. The parame- 
ters used were those typical for the experiment reported in 
Ref. 1^; uj/2tv = 70 Hz is the geometric mean of the trap 
frequencies, m is the mass of potassium-4Q, 7 — 1/350 Hz, 
and the scattering is dominated by the triplet channel with 
a — 157ao where ao is the Bohr radius ^]. The initial condi- 
tions used are given in Fig. |^. (a) The optimized cut energy 
K divided by the Fermi energy Ep. For a three-dimensional 
harmonic oscillator the Fermi energy is 

The number of atoms in each component, (c) The tempera- 
ture T, defined for the truncated Fermi-Dirac distribution as 
T = l/(/3fc_B). The temperature is normalized by the Fermi 
temperature Tf = Ep/kB- (d) The chemical potential p in 
units of Ef- (e) The Fermi temperature Tf- 

This is shown by the elastic collision rate F defined by 
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F(ep)F(e,)[l-F(e™)][l-F(e„)], 
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which is illustrated in Fig. ||. As the chemical potential 
becomes positive, and Pauli blocking of available final 
states begins to play an important role, the elastic colli- 
sion rate falls sharply. Since the elastic collision rate de- 
termines the timescale for rethermalization, at the end of 
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the simulation, evaporative cooling has virtually ceased. 
As this figure dramatically illustrates, towards the end of 
the evaporation trajectory, the elastic collision rate may 
be more than an order of magnitude suppressed from the 
value it would have if Pauli blocking of final states was 
absent. 
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FIG. 3. Elastic collision rate per atom for the optimized 
evaporation trajectory. The solid line shows the true colli- 
sion rate as defined in Eq. (|l^). The dashed line shows what 
the collision rate would be in the absence of Pauli blocking, 
by artificially replacing the [1 — _F(em)][l — F{e„)] factors in 
Eq. ^ by unity. 

In Fig. ^, we illustrate this trajectory on a semilog 
graph of temperature versus number and compare with 
experimental data. The data are taken by evaporating 
a two-component gas of as described in Ref. 
For the portion of the evaporation trajectory shown in 
Fig. ^ evaporation occurs using a 50/50 mixture of two 
spin states confined in a cylindrically symmetric har- 
monic trap whose radial frequency is 137 Hz and axial 
frequency is 19.5 Hz. After evaporation, one of the spin 
components is removed quickly (within 0.3 s) with the 
application of a frequency-swept microwave field; this re- 
moval provides a small amount of additional evaporative 
cooling that reduces the cloud temperature by 20%. The 
comparison shown in Fig. ^ illustrates that although it 
is possible theoretically to reach very low temperatures, 
the data corresponds to less efficient evaporative cool- 
ing and is presumably limited by experimental artifacts. 
These experimental limitations could include heating of 
the trapped gas, finite energy resolution of the evapora- 
tive cut, reduced dimensionality of the evaporative cut, 
and other similar problems. 
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FIG. 4. Optimum evaporation trajectory. The □ symbols 
show points on the theoretical evaporation trajectory at time 
intervals of 5 s. For comparison, the • symbols show exper- 
imental data points giving a typical evaporation trajectory. 
The A symbols show the lowest T/Tf achieved in our cur- 
rent experiment, as described in the text. 

In the current experiment we have recently found that 
improving the stability of the magnetic trapping field in- 
creased the highest achievable quantum degeneracy from 
T/Tf = 0.5 to T/Tf = 0.3. The low temperature part 
of an experimental exaporation trajectory that reached 
T/Tf = 0.3 is also shown in Fig. |^. For this data we 
used a much slower removal of the second spin compo- 
nent (within 25 s) to provide additional evaporative cool- 
ing which is not included in the theory. The experimen- 
tal progress suggests that further technical improvements 
may enable experiments to approach the low T/Tf values 
that appear possible theoretically. Furthermore Fig. ^ 
shows that the dramatic suppression of the elastic colli- 
sion rate due to Pauli blocking could be observed at the 
lowest temperatures of current experiments. 
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